function diff = moments(theta,param,aux)

    
    %% Transform parameters
    param = param_trans( theta, param);
    
    %% Solve model given parameters
    param = model( param, aux);
  
    %% Calculate inaction rates
    if  param.K>0

        T_y = KFE_firms( param, aux);               
        ln_m=linspace(log(param.m_l),log(param.m_u),aux.n_m)';
        T_yy=T_y.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        T_yy=T_yy.*repmat((exp(ln_m')>param.m_h).*(exp(ln_m')<param.m_e),aux.n_m,1);

        q   = q_fun(exp(ln_m),param);             
        q_p = q(2:end)-q(1:end-1);
        h1  = [q_p;0]./sum(q_p);
        Tr  = expm(T_yy);
        Tr  = Tr.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        inac= sum(Tr^12*h1);
        
    end
    
    %% Calculate remaining moments 
    
    u = u_fun(param);%unemployment rate     
    m = linspace(param.m_l,param.m_u,aux.n)'; 
    q = q_fun( m, param);             
    q_p = q(2:end)-q(1:end-1);
    q_p = q_p./sum(q_p);    
    job_to_job= jtj_num(delta_fun(m,param),q_fun(m,param));
    E_w = q_p'*wage_fun(m(2:end),param);% averate wage

    delta = delta_fun(m,param);
	delta_p=-(delta(2:end)-delta(1:end-1));

    mm=repmat(m(1:end-1),1,length(m(1:end-1)));
    log_gain=sum(sum((delta_p.*q_p'.* tril(log(wage_fun(mm,param)./wage_fun(mm',param))))))./sum(sum(tril((delta_p.*q_p'))));% log wage gain
    
    %% Calculate difference between moments and targets
    % Vector of differences 
    vec = [10*(job_to_job-param.job_to_job_mom),10*(u-param.u_mom), param.c./E_w-param.c_mom, log_gain - param.lnw_gain_mom];
    
    if param.K>0
       vec = [vec, inac-param.inac_mom];
    end     
    
    diff=vec*vec';
    
    
end